%matplotlib inline
import sys
sys.path.insert(0,'..')
from IPython.display import HTML,Image,SVG,YouTubeVideo
from helpers import header
HTML(header())
import matplotlib.pyplot as plt
import numpy as np
from skimage.data import imread
from skimage.morphology import erosion
def hit_or_miss(X,B12):
B1 = B12 == 1
B2 = B12 == 0
r = np.logical_and(erosion(X,B1),erosion(1-X,B2))
return r
def rotate_4(B):
#returns structuring element oriented in 4 directions
return [B,np.rot90(B),np.rot90(np.rot90(B)),np.rot90(np.rot90(np.rot90(B)))]
X = (imread('http://homepages.ulb.ac.be/~odebeir/data/man.tif')>0)[:,:,0].astype(np.uint8)
B = np.array([[2,1,2],[0,1,1],[0,0,2]]) # . are coded with 2
#iterate on the four structuring element
selem = rotate_4(B)
R = X.copy()
for i in range(10):
for s in selem:
HoM = hit_or_miss(R,s)
R[HoM] = 0
plt.figure(figsize=[10,10])
plt.subplot(1,2,1)
plt.imshow(X,interpolation='nearest')
plt.subplot(1,2,2)
plt.imshow(X,interpolation='nearest',alpha = .8)
plt.imshow(R,interpolation='nearest',alpha=.5);
A special case of thinning is the skeleton extraction, by applying the thickening for these two set of structuring elements, one can reduce every $X$ to single line structures sharing the same topology.
$$ \begin{matrix} 0 & 0 & 0 \\ \cdot & 1 & \cdot \\ 1 & 1 & 1 \end{matrix} $$and
$$ \begin{matrix} \cdot & 0 & 0 \\ 1 & 1 & 0 \\ \cdot & 1 & \cdot \end{matrix} $$
and there 3 other orientations.
# skeletoninzing structuring elements
B1 = np.array([[0,0,0],[2,1,2],[1,1,1]])
B2 = np.array([[2,0,0],[1,1,0],[2,1,2]])
selem = rotate_4(B1) + rotate_4(B2) #join the two list of structuring element
X[1,15:17] = 1
R = X.copy()
for i in range(4):
for s in selem:
HoM = hit_or_miss(R,s)
R[HoM] = 0
plt.figure(figsize=[10,10])
plt.subplot(1,2,1)
plt.imshow(X,interpolation='nearest')
plt.subplot(1,2,2)
plt.imshow(X,interpolation='nearest',alpha = .8)
plt.imshow(R,interpolation='nearest',alpha=.5);
an application of the skeleton extraction is found in the pre-processing of fingerprints where we are looking for specific minutiae
from skimage.filter import threshold_otsu
im = imread('http://biometrics.idealtest.org/userfiles/image/FingerprintV5Fig.1.jpg')[:,:150,0]
th_otsu = threshold_otsu(im)
R = (im<th_otsu).copy()
for i in range(4):
for s in selem:
HoM = hit_or_miss(R,s)
R[HoM] = 0
plt.figure(figsize=[10,5])
plt.subplot(1,2,1)
plt.imshow(im,cmap=plt.cm.gray);
plt.subplot(1,2,2)
plt.imshow(R,cmap=plt.cm.gray);
/home/olivier/anaconda/lib/python2.7/site-packages/skimage/filter/__init__.py:6: skimage_deprecation: The `skimage.filter` module has been renamed to `skimage.filters`. This placeholder module will be removed in v0.13. warn(skimage_deprecation('The `skimage.filter` module has been renamed '
see also:
The skeleton, or other morphological processing can produce unwanted small artifacts, these irragularities can often be removed used by a process called pruning.
Pruning is a specific thinning algorithm using the following structuring element (and ther rotations):
$$ \begin{matrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & \cdot & \cdot \end{matrix} $$and
$$ \begin{matrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ \cdot & \cdot & 0 \end{matrix} $$
and there 3 other orientations, example of a pruning applied iteratively on a binary skeleton:
ima = imread('http://homepages.ulb.ac.be/~odebeir/data/programs.png')[:,:,0] >10
# skeletoninzing structuring elements
B1 = np.array([[0,0,0],[2,1,2],[1,1,1]])
B2 = np.array([[2,0,0],[1,1,0],[2,1,2]])
selem = rotate_4(B1) + rotate_4(B2) #join the two list of structuring element
R = ima.copy()
for i in range(12):
for s in selem:
HoM = hit_or_miss(R,s)
R[HoM] = 0
plt.figure(figsize=[10,10])
plt.subplot(1,2,1)
plt.imshow(ima,interpolation='nearest',cmap=plt.cm.gray)
plt.subplot(1,2,2)
plt.imshow(ima,interpolation='nearest',alpha = .8,cmap=plt.cm.gray)
plt.imshow(R,interpolation='nearest',alpha=.5,cmap=plt.cm.gray);
# pruning structuring elements
B1 = np.array([[0,0,0],[0,1,0],[0,2,2]])
B2 = np.array([[0,0,0],[0,1,0],[2,2,0]])
selem = rotate_4(B1) + rotate_4(B2) #join the two list of structuring element
plt.figure(figsize=[15,25])
k = 1
for i in range(24):
for s in selem:
HoM = hit_or_miss(R,s)
R[HoM] = 0
if k%3 == 0:
plt.subplot(4,2,k/3)
plt.imshow(ima,interpolation='nearest',alpha = .7,cmap=plt.cm.gray)
plt.imshow(R.copy(),interpolation='nearest',alpha=.8,cmap=plt.cm.gray);
k += 1
The distance function $dist_X$ associate to a set $X$, the distance for each point of $X$ to the background
$$ dist_X(p) = min\{n \in \mathbb \; | \; p \notin X \ominus nB \}$$explicit distance map (brute force) computing can be computer intensive, therefore other appriximated distance are proposed
Image('http://homepages.ulb.ac.be/~odebeir/data/distances.png')
from scipy.ndimage import distance_transform_cdt,distance_transform_edt,distance_transform_bf
import matplotlib.pyplot as plt
import numpy as np
X = np.asarray([[0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0],
[0,0,0,1,1,1,0,0,0],
[0,0,1,1,1,0,0,0,0],
[0,0,1,1,1,1,0,0,0],
[0,0,1,1,1,0,0,0,0],
[0,0,0,1,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0]])
dist_X = distance_transform_edt(X)
dist_bf = distance_transform_bf(X)
dist_c = distance_transform_cdt(X)
plt.figure(figsize=[8,8])
plt.subplot(2,2,1)
plt.imshow(X,interpolation='nearest',cmap=plt.cm.gray);X = np.zeros((10,10))
plt.title('X')
plt.colorbar()
plt.subplot(2,2,2)
plt.imshow(dist_X,interpolation='nearest',cmap=plt.cm.gray);X = np.zeros((10,10))
plt.title('Euclidian distance map')
plt.colorbar();
plt.subplot(2,2,3)
plt.imshow(dist_c,interpolation='nearest',cmap=plt.cm.gray);X = np.zeros((10,10))
plt.title('Chamfer distance map')
plt.colorbar()
plt.subplot(2,2,4)
plt.imshow(dist_bf,interpolation='nearest',cmap=plt.cm.gray);X = np.zeros((10,10))
plt.title('Brute force distance map')
plt.colorbar();
import numpy as np
from skimage.morphology import disk
import skimage.filter.rank as skr
from scipy import ndimage as ndi
np.random.seed(1)
n = np.random.random((256,256))<.0005
d = skr.maximum(n,disk(15))
distance = ndi.distance_transform_edt(d).astype(np.uint8)
local_max = distance == skr.maximum(distance,disk(20))
background = distance == 0
c_local_max = local_max.copy()
c_local_max[background] = 0
plt.figure(figsize=[13,8])
plt.subplot(2,3,1)
plt.imshow(d)
plt.subplot(2,3,2)
plt.imshow(distance)
plt.subplot(2,3,3)
plt.imshow(skr.maximum(local_max,disk(2)),cmap=plt.cm.gray)
plt.subplot(2,3,4)
plt.imshow(background)
plt.subplot(2,3,5)
plt.imshow(d,alpha=.5)
plt.imshow(skr.maximum(c_local_max,disk(2)),cmap=plt.cm.gray,alpha=.5);
see also:
The object labeling consist in assigning a unique integer value to each connected component of a binary image.
from skimage.measure import label
ima = imread('http://homepages.ulb.ac.be/~odebeir/data/programs.png')[:,:,0] >10
lab = label(ima,background=0)
plt.figure(figsize=[20,5])
plt.subplot(1,2,1)
plt.imshow(ima)
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(lab)
plt.colorbar();